{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Using crack submodels in PyBaMM\n",
    "In this notebook we show how to use the crack submodel with battery DFN or SPM models. To see all of the models and submodels available in PyBaMM, please take a look at the documentation [here](https://docs.pybamm.org/en/latest/source/api/models/index.html)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-09-26T23:15:29.863147Z",
     "start_time": "2023-09-26T23:15:29.848113Z"
    }
   },
   "outputs": [],
   "source": [
    "%pip install \"pybamm[plot,cite]\" -q    # install PyBaMM if it is not installed\n",
    "import os\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import pybamm\n",
    "\n",
    "os.chdir(pybamm.__path__[0] + \"/..\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then we load the DFN, SPMe or SPM, by choosing one and commenting the others. \n",
    "\n",
    "When you load a model in PyBaMM it builds by default. Building the model sets all of the model variables and sets up any variables which are coupled between different submodels: this is the process which couples the submodels together and allows one submodel to access variables from another. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-09-26T23:15:29.896170Z",
     "start_time": "2023-09-26T23:15:29.885253Z"
    }
   },
   "outputs": [],
   "source": [
    "model = pybamm.lithium_ion.DFN(\n",
    "    options={\n",
    "        \"particle\": \"Fickian diffusion\",\n",
    "        \"particle mechanics\": \"swelling and cracking\",  # other options are \"none\", \"swelling only\"\n",
    "    }\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Load the parameter set Ai2020 which contains mechanical parameters. Other sets may not contain mechanical parameters should add them manually. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-09-26T23:15:29.896383Z",
     "start_time": "2023-09-26T23:15:29.887700Z"
    }
   },
   "outputs": [],
   "source": [
    "param = pybamm.ParameterValues(\"Ai2020\")\n",
    "## It can update the speed of crack propagation using the commands below:\n",
    "# param.update({\"Negative electrode Cracking rate\":3.9e-20*10}, check_already_exists=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can get the default parameters for the model and update them with the parameters required by the cracking model. Eventually, we would like these to be added to their own chemistry (you might need to adjust the path to the parameters file to your system).\n",
    "Now the model can be processed and solved in the usual way, and we still have access to model defaults such as the default geometry and default spatial methods"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "Depending on the parameter set being used, the particle cracking model can require a large number of mesh points inside the particles to be numerically stable."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-09-26T23:15:29.901590Z",
     "start_time": "2023-09-26T23:15:29.893646Z"
    }
   },
   "outputs": [],
   "source": [
    "var_pts = {\n",
    "    \"x_n\": 20,  # negative electrode\n",
    "    \"x_s\": 20,  # separator\n",
    "    \"x_p\": 20,  # positive electrode\n",
    "    \"r_n\": 26,  # negative particle\n",
    "    \"r_p\": 26,  # positive particle\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-09-26T23:15:31.412965Z",
     "start_time": "2023-09-26T23:15:30.060671Z"
    }
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "70c1339268f44ef19ee852fe5bd01653",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(FloatSlider(value=0.0, description='t', max=1.0, step=0.01), Output()), _dom_classes=('w…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "sim = pybamm.Simulation(\n",
    "    model,\n",
    "    parameter_values=param,\n",
    "    var_pts=var_pts,\n",
    ")\n",
    "solution = sim.solve(t_eval=[0, 3600], inputs={\"C-rate\": 1})\n",
    "# plot\n",
    "quick_plot = pybamm.QuickPlot(solution)\n",
    "quick_plot.dynamic_plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the results as required."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-09-26T23:15:31.665271Z",
     "start_time": "2023-09-26T23:15:31.432275Z"
    }
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "51b97221822843719c6bea4506288ed3",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(FloatSlider(value=0.0, description='t', max=3600.0, step=10.0), Output()), _dom_classes=…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# extract voltage\n",
    "E_n = param[\"Negative electrode Young's modulus [Pa]\"]\n",
    "stress_t_n_surf = solution[\"Negative particle surface tangential stress [Pa]\"]\n",
    "x = solution[\"x [m]\"].entries[0:19, 0]\n",
    "c_s_n = solution[\"Negative particle concentration\"]\n",
    "r_n = solution[\"r_n [m]\"].entries[:, 0, 0]\n",
    "\n",
    "# plot\n",
    "\n",
    "\n",
    "def plot_concentrations(t):\n",
    "    f, (ax1, ax2, ax3, ax4) = plt.subplots(1, 4, figsize=(20, 4))\n",
    "    ax1.plot(x, stress_t_n_surf(t=t, x=x) / E_n)\n",
    "    ax1.set_xlabel(r\"$x_n$ [m]\")\n",
    "    ax1.set_ylabel(\"$\\sigma_t/E_n$\")\n",
    "\n",
    "    (plot_c_n,) = ax2.plot(\n",
    "        r_n, c_s_n(r=r_n, t=t, x=x[0])\n",
    "    )  # can evaluate at arbitrary x (single representative particle)\n",
    "    ax2.set_ylabel(\"Negative particle concentration\")\n",
    "    ax2.set_xlabel(r\"$r_n$ [m]\")\n",
    "    ax2.set_ylim(0, 1)\n",
    "    ax2.set_title(\"Close to current collector\")\n",
    "    ax2.grid()\n",
    "\n",
    "    (plot_c_n,) = ax3.plot(\n",
    "        r_n, c_s_n(r=r_n, t=t, x=x[10])\n",
    "    )  # can evaluate at arbitrary x (single representative particle)\n",
    "    ax3.set_ylabel(\"Negative particle concentration\")\n",
    "    ax3.set_xlabel(r\"$r_n$ [m]\")\n",
    "    ax3.set_ylim(0, 1)\n",
    "    ax3.set_title(\"In the middle\")\n",
    "    ax3.grid()\n",
    "\n",
    "    (plot_c_n,) = ax4.plot(\n",
    "        r_n, c_s_n(r=r_n, t=t, x=x[-1])\n",
    "    )  # can evaluate at arbitrary x (single representative particle)\n",
    "    ax4.set_ylabel(\"Negative particle concentration\")\n",
    "    ax4.set_xlabel(r\"$r_n$ [m]\")\n",
    "    ax4.set_ylim(0, 1)\n",
    "    ax4.set_title(\"Close to separator\")\n",
    "    ax4.grid()\n",
    "    plt.show()\n",
    "\n",
    "\n",
    "import ipywidgets as widgets\n",
    "\n",
    "widgets.interact(\n",
    "    plot_concentrations, t=widgets.FloatSlider(min=0, max=3600, step=10, value=0)\n",
    ");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot results using the default functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-09-26T23:15:31.932322Z",
     "start_time": "2023-09-26T23:15:31.685338Z"
    }
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "3a5afe2d4d7a4987af6c97d9f57b5872",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(FloatSlider(value=0.0, description='t', max=1.0, step=0.01), Output()), _dom_classes=('w…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "label = [\"Crack model\"]\n",
    "output_variables = [\n",
    "    \"Negative particle crack length [m]\",\n",
    "    \"Positive particle crack length [m]\",\n",
    "    \"X-averaged negative particle crack length [m]\",\n",
    "    \"X-averaged positive particle crack length [m]\",\n",
    "]\n",
    "quick_plot = pybamm.QuickPlot(\n",
    "    solution, output_variables, label, variable_limits=\"tight\"\n",
    ")\n",
    "quick_plot.dynamic_plot();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## References\n",
    "\n",
    "The relevant papers for this notebook are:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-09-26T23:15:31.959902Z",
     "start_time": "2023-09-26T23:15:31.946562Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1] Weilong Ai, Ludwig Kraft, Johannes Sturm, Andreas Jossen, and Billy Wu. Electrochemical thermal-mechanical modelling of stress inhomogeneity in lithium-ion pouch cells. Journal of The Electrochemical Society, 167(1):013512, 2019. doi:10.1149/2.0122001JES.\n",
      "[2] Joel A. E. Andersson, Joris Gillis, Greg Horn, James B. Rawlings, and Moritz Diehl. CasADi – A software framework for nonlinear optimization and optimal control. Mathematical Programming Computation, 11(1):1–36, 2019. doi:10.1007/s12532-018-0139-4.\n",
      "[3] Rutooj Deshpande, Mark Verbrugge, Yang-Tse Cheng, John Wang, and Ping Liu. Battery cycle life prediction with coupled chemical degradation and fatigue mechanics. Journal of the Electrochemical Society, 159(10):A1730, 2012. doi:10.1149/2.049210jes.\n",
      "[4] Marc Doyle, Thomas F. Fuller, and John Newman. Modeling of galvanostatic charge and discharge of the lithium/polymer/insertion cell. Journal of the Electrochemical society, 140(6):1526–1533, 1993. doi:10.1149/1.2221597.\n",
      "[5] Charles R. Harris, K. Jarrod Millman, Stéfan J. van der Walt, Ralf Gommers, Pauli Virtanen, David Cournapeau, Eric Wieser, Julian Taylor, Sebastian Berg, Nathaniel J. Smith, and others. Array programming with NumPy. Nature, 585(7825):357–362, 2020. doi:10.1038/s41586-020-2649-2.\n",
      "[6] Valentin Sulzer, Scott G. Marquis, Robert Timms, Martin Robinson, and S. Jon Chapman. Python Battery Mathematical Modelling (PyBaMM). Journal of Open Research Software, 9(1):14, 2021. doi:10.5334/jors.309.\n"
     ]
    }
   ],
   "source": [
    "pybamm.print_citations()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.18"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": true
  },
  "vscode": {
   "interpreter": {
    "hash": "187972e187ab8dfbecfab9e8e194ae6d08262b2d51a54fa40644e3ddb6b5f74c"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
